\(Y_t\) is ARMA with error \(\epsilon_t\). \(\epsilon_t\) is GARCH. \[ \begin{align} \Phi(B) Y_t &= \Theta(B) \epsilon_t \\ \\ \epsilon_t &= \sigma_t e_t \hspace10mm e_t \sim_{iid} N(0,1) \\\\ \sigma_t^2 &= \omega + \alpha \epsilon_{t-1}^2 + \beta \sigma^2_{t-1} \end{align} \]
Daily Price of SP500 ETF (SPY) from Jan 02 2000 to Dec 31 2014
## Warning in if (!header) rlabp <- FALSE: the condition has length > 1 and
## only the first element will be used
## Warning in if (header) {: the condition has length > 1 and only the first
## element will be used
## ixDay Date Weekday X SPY.Open SPY.High SPY.Low SPY.Close
## 1 1 1/3/2000 Monday 1 148.25 148.25 143.88 145.44
## 2 2 1/4/2000 Tuesday 2 143.53 144.06 139.64 139.75
## 3 3 1/5/2000 Wednesday 3 139.94 141.53 137.25 140.00
## 4 4 1/6/2000 Thursday 4 139.62 141.50 137.75 137.75
## 5 5 1/7/2000 Friday 5 140.31 145.75 140.06 145.75
## 6 6 1/10/2000 Monday 6 146.25 146.91 145.03 146.25
## SPY.Volume SPY.Adjusted
## 1 8164300 110.33
## 2 8089800 106.01
## 3 12177900 106.20
## 4 6227200 104.50
## 5 8066500 110.57
## 6 5741700 110.94
layout(1)
# Fit GARCH model
library(fGarch)
Fit1 <- garchFit(~ garch(1,1), data=X1, cond.dist="norm", include.mean = FALSE, trace = FALSE)
## Fit1@fit$par #-- estimated parameters
## Fit1@residuals #-- this is not the garch residuals!!!! same as X1
## Fit1@sigma.t #-- estimated sig_t
## Fit1@residuals/Fit1@sigma.t #-- this is the (standardized) GARCH residuals
print(Fit1@fit$ics) #-- AIC and BIC are here## AIC BIC SIC HQIC
## -6.470319 -6.465360 -6.470321 -6.468556
plot(X1)
lines( 1.96*Fit1@sigma.t, type="l", col="red")
lines(-1.96*Fit1@sigma.t, type="l", col="red")## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series
## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0 0 0 0.719 0.855 0 1.001
library(xts)
#--- ARMA-GARCH estimation example ---
D <- read.csv("http://gozips.uakron.edu/~nmimoto/pages/datasets/SPY.csv", header=T)## Warning in if (!header) rlabp <- FALSE: the condition has length > 1 and
## only the first element will be used
## Warning in if (header) {: the condition has length > 1 and only the first
## element will be used
## ixDay Date Weekday X SPY.Open SPY.High SPY.Low SPY.Close
## 1 1 1/3/2000 Monday 1 148.25 148.25 143.88 145.44
## 2 2 1/4/2000 Tuesday 2 143.53 144.06 139.64 139.75
## 3 3 1/5/2000 Wednesday 3 139.94 141.53 137.25 140.00
## 4 4 1/6/2000 Thursday 4 139.62 141.50 137.75 137.75
## 5 5 1/7/2000 Friday 5 140.31 145.75 140.06 145.75
## 6 6 1/10/2000 Monday 6 146.25 146.91 145.03 146.25
## SPY.Volume SPY.Adjusted
## 1 8164300 110.33
## 2 8089800 106.01
## 3 12177900 106.20
## 4 6227200 104.50
## 5 8066500 110.57
## 6 5741700 110.94
## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series
## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0 0 0 0 0 0 0.012
library(fGarch)
#- Estimae parameters of ARMA and GARCH at the same time
Fit1 <- garchFit(~ arma(1,1) + garch(1,1), data=X1, cond.dist="norm", include.mean = FALSE, trace = FALSE)
summary(Fit1) ##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~arma(1, 1) + garch(1, 1), data = X1, cond.dist = "norm",
## include.mean = FALSE, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ arma(1, 1) + garch(1, 1)
## <environment: 0x7f8ff5b69390>
## [data = X1]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## ar1 ma1 omega alpha1 beta1
## -1.4401e-01 2.4903e-01 1.5994e-06 1.0580e-01 8.8436e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## ar1 -1.440e-01 1.317e-01 -1.093 0.2743
## ma1 2.490e-01 1.302e-01 1.913 0.0557 .
## omega 1.599e-06 2.683e-07 5.961 2.51e-09 ***
## alpha1 1.058e-01 9.451e-03 11.195 < 2e-16 ***
## beta1 8.844e-01 9.345e-03 94.635 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 12225.18 normalized: 3.241034
##
## Description:
## Wed Apr 11 14:54:21 2018 by user:
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 2177.437 0
## Shapiro-Wilk Test R W 0.9765013 0
## Ljung-Box Test R Q(10) 11.66489 0.3081103
## Ljung-Box Test R Q(15) 16.35115 0.3591053
## Ljung-Box Test R Q(20) 23.32912 0.2729244
## Ljung-Box Test R^2 Q(10) 8.715104 0.5593378
## Ljung-Box Test R^2 Q(15) 11.82243 0.692417
## Ljung-Box Test R^2 Q(20) 13.78888 0.8410369
## LM Arch Test R TR^2 8.749246 0.7241858
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -6.479417 -6.471151 -6.479420 -6.476478
##Fit1@fit$par #-- estimated parameters
##Fit1@residuals #-- this is ARCH residuals. (before GARCH)
##print(Fit1@fit$ics) #-- AIC and BIC are here
##Randomness.tests(res1)
##plot(Fit1) #-- choose option 13 for residual qq plot
Fit2 <- garchFit(~ arma(0,1) + garch(1,1), data=X1, cond.dist="norm", include.mean = FALSE, trace = FALSE)
summary(Fit2) ##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~arma(0, 1) + garch(1, 1), data = X1, cond.dist = "norm",
## include.mean = FALSE, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ arma(0, 1) + garch(1, 1)
## <environment: 0x7f8ff8faae68>
## [data = X1]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## ma1 omega alpha1 beta1
## 1.0525e-01 1.5932e-06 1.0567e-01 8.8453e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## ma1 1.053e-01 1.808e-02 5.821 5.84e-09 ***
## omega 1.593e-06 2.676e-07 5.953 2.63e-09 ***
## alpha1 1.057e-01 9.442e-03 11.191 < 2e-16 ***
## beta1 8.845e-01 9.333e-03 94.771 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 12224.63 normalized: 3.240888
##
## Description:
## Wed Apr 11 14:54:22 2018 by user:
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 2127.103 0
## Shapiro-Wilk Test R W 0.9767988 0
## Ljung-Box Test R Q(10) 14.02298 0.1719459
## Ljung-Box Test R Q(15) 18.63737 0.2306394
## Ljung-Box Test R Q(20) 25.79804 0.1725862
## Ljung-Box Test R^2 Q(10) 8.686103 0.56213
## Ljung-Box Test R^2 Q(15) 11.83827 0.6912271
## Ljung-Box Test R^2 Q(20) 13.80512 0.8402378
## LM Arch Test R TR^2 8.718929 0.7267317
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -6.479655 -6.473043 -6.479657 -6.477304
Fit3 <- garchFit(~ arma(0,2) + garch(1,1), data=X1, cond.dist="norm", include.mean = FALSE, trace = FALSE)
summary(Fit3) ##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~arma(0, 2) + garch(1, 1), data = X1, cond.dist = "norm",
## include.mean = FALSE, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ arma(0, 2) + garch(1, 1)
## <environment: 0x7f8ff2f6c548>
## [data = X1]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## ma1 ma2 omega alpha1 beta1
## 1.0242e-01 -2.7523e-02 1.5971e-06 1.0548e-01 8.8469e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## ma1 1.024e-01 1.771e-02 5.782 7.37e-09 ***
## ma2 -2.752e-02 1.761e-02 -1.563 0.118
## omega 1.597e-06 2.681e-07 5.958 2.56e-09 ***
## alpha1 1.055e-01 9.429e-03 11.187 < 2e-16 ***
## beta1 8.847e-01 9.326e-03 94.861 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 12226.29 normalized: 3.241329
##
## Description:
## Wed Apr 11 14:54:22 2018 by user:
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 2228.007 0
## Shapiro-Wilk Test R W 0.9761852 0
## Ljung-Box Test R Q(10) 10.67642 0.3832771
## Ljung-Box Test R Q(15) 15.33805 0.4273523
## Ljung-Box Test R Q(20) 21.88992 0.3465134
## Ljung-Box Test R^2 Q(10) 8.906359 0.541013
## Ljung-Box Test R^2 Q(15) 11.88314 0.6878524
## Ljung-Box Test R^2 Q(20) 13.86514 0.8372674
## LM Arch Test R TR^2 8.853235 0.7154123
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -6.480008 -6.471742 -6.480011 -6.477069
sig.t <- xts(Fit2@sigma.t, order.by=index(X1)) #-- estimated sig_t fo GARCH
res1 <- xts(Fit2@residuals/Fit2@sigma.t, order.by=index(X1)) #-- this is the (standardized) ARMA-GARCH residuals
Randomness.tests(as.numeric(res1))## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series
## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0.231 0.173 0.172 0.691 0.84 0 1.001
#- plot log-return and sigma_t -
layout(matrix(1:2, 2, 1, byrow=T))
plot(X1, main="SPY daily log-return")
plot(sig.t, type="l", main="estimated sigma_t") #- estimated Sigma_t #--- h-step ahead prediction from ARMA-GARCH ---
h = 15 #- predict 15 days ahead
X.pred <- xts(predict(Fit2, n.ahead=h)[,1], order.by=index(X1)[length(X1)]+(1:h))
SD.pred <- xts(predict(Fit2, n.ahead=h)[,3], order.by=index(X1)[length(X1)]+(1:h))
#- plot log-return overlay with sigma_t -
plot(rbind(X1["2014"], X.pred))
lines(X.pred, col="red")
lines( 1.96*SD.pred, type="l", col="red")
lines(-1.96*SD.pred, type="l", col="red") #--- Rolling 1-step prediction from ARMA-GARCH
Y <- X1
Y.pred1 <- sig.pred1 <- date.stp <- numeric(0)
for (i in 1:1000){
T <- Y[(1:250)+(i-1)]
Fit44 <- garchFit(~ arma(0,1) + garch(1,1), data=T, cond.dist="norm", include.mean = FALSE, trace = FALSE)
Y.pred1[i] <- predict(Fit44, n.ahead=1)[,1]
sig.pred1[i] <- predict(Fit44, n.ahead=1)[,3]
date.stp[i] <- index(Y[250+i])
}## Warning in sqrt(diag(fit$cvar)): NaNs produced
Y.pred <- xts(Y.pred1, order.by=as.Date(date.stp))
sig.pred <- xts(sig.pred1, order.by=as.Date(date.stp))
#- Compare predicted vs actual -
plot(Y["2000::2004"]); lines(Y.pred, col="red"); lines( 1.96*sig.pred, col="red"); lines(-1.96*sig.pred, col="red")## ..1 ..2
## 2000-12-29 -0.0106150060 0.0009925981
## 2001-01-02 -0.0280598120 -0.0012141588
## 2001-01-03 0.0007836377 -0.0029252764
## 2001-01-04 0.0409786484 0.0003025638
## 2001-01-05 -0.0290649392 0.0040535907
## 2001-01-08 -0.0116787375 -0.0026158966
## [1] 0.53
## [1] 0.508
## [1] 0.33979142 -0.28412857 -0.15940948 -0.07124178 0.43093553
## [6] -0.30531458 1.89158391 1.57043769 1.26749188 -0.38981266
In iid setting, K-S test can be used to test hypothesis \(X_i \sim F\), and
\[ H_0: F = F_0 \hspace5mm vs \hspace5mm H_A: F \ne F_0 \hspace10mm \mbox{{\scriptsize ($F_0$ is completely specified) } } \]
K-S test uses Empirical Distibution (EDF, ECDF), and calculate the test statistic \[ KS \hspace3mm = \hspace3mm \sup_x |\hat F_n(x) - F(x) | \]
Where is \((\frac{i}{n}, X_{(i)})\) and \((\frac{i-1}{n}, X_{(i)})\)?
\[ \begin{align} D_n &= \sup_x |\hat F_n(x) - F_0(x)| \\\\ &= \max_i \Big( \Big|\frac{i}{n} - F_0(x_{(i)}) \Big| , \hspace3mm \Big|\frac{i-1}{n} - F_0(x_{(i)}) \Big| \Big) \\ \\ \end{align} \]
If \(X_i\) are indeed from \(F_0\), \[ D_n \hspace3mm =_d \hspace3mm \max_i \Big( \Big|\frac{i}{n} - U_{(i)} \Big| , \hspace3mm \Big|\frac{i-1}{n} - U_{(i)} \Big| \Big) \] with \(U_i \sim_{iid} Unif(0,1)\).
No matter what \(F\) is, distribution of the test statistic KS under the null is known.
KS test with completely specified \(F_0\) is a distribution-free test.
When \(F_0\) contains nuisance parameter, Distribution of \(D_n\) under the null depends on \(F\). KS test becomes only asymptotically distribution free test. For finite \(n\), null distribution must be computed by Monte Carlo simulation.
In time sries, we want to use K-S test for checking distribution of innovations. (errors).
Since we assume \(e_t\) are i.i.d. noise, if we can observe \(e_t\), K-S test can be directly applicable.
However, in Time Series analysis \(e_t\) are not observable, and only residuals \(\hat e_t\) are available.
Can we still use K-S test on \(\hat e_t\)?